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Abstract 

Molecular Dynamics simulations were carried out to test the ther- 
modynamic theory of weakly excited, two-dimensional, granular sys- 
tems [Hayakawa and Hong, Phys. Rev. Lett. 78, 2764(1997)], where 
granular materials are viewed as a collection of spinless Fermions. We 
first determine the global temperature T by fitting the steady state 
density profile to the Fermi distribution function, and then measure 
the center of mass, < z(T) >, and its fluctuations, < (Az(T)) 2 > as 
a function of T. We find a fairly good agreement between theory and 
simulations, in particular, in the estimation of the temperature and 
the scaling behavior of < z(T) > and < (Az(T)) 2 >. 

PACS numbers: 81.05 Rm, 05., 07.05Tp, 82.20.wt 
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I. Introduction 

In a recent paper, Hayakawa and Hong(HH) [1] explored a simple con- 
sequence of the excluded volume effect for a dense granular system. They 
proposed a thermodynamic theory for weakly excited, dissipative, nonequi- 
librium, granular systems from the view point of the elementary excitation 
of non-interacting spinless Fermions. Based on a simple thermodynamic 
argument and combined with the observation that the granular state in a 
vibrating bed may be viewed as an excited state, HH demonstrated that the 
steady state density profile must be given by the Fermi distribution function. 
This enables one to associate and define a macroscopic, global parameter, T, 
which is suitable to characterize the state of the vibrating system. The pa- 
rameter T is the thermodynamic temperature, and its relation to external 
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control parameters has been derived in ref.[l]. The theory of Hayakawa and 
Hong has been successful in explaining the observed density profile of the ex- 
perimental data of Clement and Rajchenbach [2]. The purpose of this paper 
is to further test the predictions of HH using Molecular Dynamics simula- 
tions, in particular, with regards to the configurational statistics of granular 
materials in a vibrating bed. We first briefly explain the background of the 
thermodynamic theory in section II along with a summary of the theory of 
ref.[l] in section III. We then present the simulation results in section IV. 
In section V, is a discussion of some of the important issues related to this 
work. 

II. Background of Fermi Statistics and Thermodynamic Theory of 
Granular Materials 

The system being studied here is a dense, dissipative, nonequilibrium, 
granular system, where the mean free path of the grains is of the order of a few 
particle diameters. Hence, each particle may be considered to be effectively 
confined in a cage as in the free volume theory of a dense liquid [4]. In such 
a case, an observation has been made in ref.[l] that the basic granular state 
is not a gas, but a solid or crystal, and thus the effective thermodynamic 
theory based on the free energy argument may be more appropriate than 
the kinetic theory in studying this state. In such the configurational 

statistics of the steady state may be determined by the variational method 
as the most probable or minimum free energy state. 

To be more specific, consider the excitation of disordered granular mate- 
rials confined in a box with vibrations of the bottom plate. The vibrations 
will inject energy into the system which cause the ground state to become 
unstable, and a newly excited state will emerge with an expanded volume. 
The time averaged configurational states of this new excited state have un- 
dergone structural distortions. However, the degree of distortions from the 
ground state may be small for a weakly excited state, possibly justifying the 
use of an effective thermodynamic theory based on the variational principle. 
Such a thermodynamic approach may be further justified by the following 
two recent experiments: 

1. Weakly or moderately excited regime: Clement and Rajchenbach(CR) 
[2] have performed an experiment with the vibrational strength, T, of the 
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order one for a two dimensional vibrating bed, using inclined side walls to 
suppress convections. Here, T = Au 2 /g with A and u, the amplitude and 
frequency of the vibrating plane, and g, the gravitational constant. CR have 
found that the ensemble-averaged density profile as a function of height from 
the bottom layer obeys a universal function that is independent of the phase 
of oscillations of the vibrating plate. Namely, it is independent of the kinetics 
imposed on the system. One conceptually important point here is that the 
reference point of the density profile is not the bottom plate, but the bottom 
layer, which of course is fluidized. 

2. Highly excited regime: Warr and Hansen(WH) [5] have performed 
an experiment on highly agitated, vertically vibrating beds of T « 30 — 50 
using steel balls with a small coefficient of restitution. They have found 
that the collective behaviors of this vibrated granular medium in a station- 
ary nonequilibrium state exhibits strong similarities to those of an atomistic 
fluid in thermal equilibrium at the corresponding particle packing fraction, 
in particular, in the two-point correlation function [6]. 

The results of both experiments indicate that for both moderate or highly 
excited systems, one-to-one correspondence seems to exist between configu- 
rational statistics of the nonequilibrium stationary state and the equilibrium 
thermal state. In fact, this is not so surprising considering that upon vibra- 
tion, the granular materials expand and increase the volume of the system. 
In turn, this increase corresponds to a rise in the potential energy after the 
conflgurational average is appropriately taken. Then the problem reduces 
to the packing problem, and the temperature-like variable, T, can be as- 
sociated to the vibrating bed. The existence of distinctive conflgurational 
statistics in the density profile of CR (and also in WH in a special case) 
appears to be fairly convincing evidence that kinetic aspects of the excited 
granular materials may be separated out from the statistical configurations. 
These observations are the basis of the thermodynamic theory proposed in 
[1]. Note that the Fermi statistics is essentially the macroscopic manifesta- 
tion of the classical excluded volume effect and the anisotropy which causes 
the ordering of potential energy by gravity. The top surface of the granules 
plays the role of a Fermi surface, and the thin boundary layers that appear 
near the top layer upon excitation play the role of excited electrons of the 
Fermi gas in metals. 
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III. Thermodynamic Theory of Weakly Excited Granular Systems 



1. Fermi temperature: In ref.[l], the vibrating system was viewed from 
two different points of view: one may view it as a mechanical system, in which 
case the expansion is due to excitation induced by mechanical vibrations 
of vibration strength T = Au 2 /g. In this case, the expansion is purely 
due to kinetics. In an attempt to develop a thermodynamic theory, such 
a system was also viewed in ref.[l] as a thermal system in contact with a 
heat reservoir, and a global temperature T was associated with it. In this 
case, the expansion is purely a thermal expansion. By equating the thermal 
expansion, Ah, defined as the increase in the center of mass, and the kinetic 
expansion, g/u 2 H (T), where H (T) is the jump height of a single ball [1,7] 
at the vibrating plate(see eq.(l) in ref.[l]), a closure in the thermodynamic 
theory of powders was obtained in ref.[l]. Since the density decrease above 
the Fermi surface is not sharp, but smooth, one may replace gH {T)/uo 2 with 
h (T)/a, where h (T) is the maximum jump height of a single ball at the 
Fermi surface(or vibrating plate) determined by MD simulations. The factor 
a was introduced to incorporate (i) the smooth decrease in the density profile 
near the Fermi surface, and (ii) the suppression in the jump height due to 
dissipation. By equating the kinetic expansion and the thermal expansion, 
one obtains the follwing explicit relationships between the temperature T 
and the control parameters: 



Note that when a single particle is on a vibrating plate, the energy from the 
bottom wall is transferred to the particle through direct contact. In the case 
of many particles, the supplied energy at the vibrating plate must first travel 
through other particles locked in their respective lattice states before reaching 
those particles in the fluidized layer. For (la), gH (T)/uj 2 , the maximum of 
the single ball jump height, is a lower bound because the relative velocity 
between the ball and the plate is assumed to be zero, making this temperature 
a lower limit for the system. In reality the relative velocity can be much 
higher and a could be smaller than 1. For (lb), the expansion of the fluidized 
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layer near the Fermi surface is certainly less than the simulated single ball 
jump height because of the dissipation of energies through collisions. In this 
case, a is like a dissipation constant. The best fit a's used in this work were 
a = 1 in (la) and a = 64/5 (fractional fit) in (lb). The second value of a fit 
all the data, regardless of whether the sine wave or the triangular wave was 
used to vibrate the bottom wall. 



2. The center of mass and its fluctuations: Since the density profile 
is given by the Fermi function, it is straightfoward to compute the center of 
mass, z(T), and its fluctuations, < (Az) 2 >. The following formulas were 
derived in ref.[l]: 

MT) = ,<r> - ,(o) = ^fi^f (2«) 



< (Azf >=< (z(T)- < z >f >=< (Ah) 2 > I A = -(—f—. (26) 

Note that the total expansion, Ah(T) = fi Q Az and its fluctuations < (Ah) 2 > 
j D 2 =< fi (Az) 2 / D 2 > are only a function of the dimensionless Fermi tem- 
perature Tf = T/mgD as expected. Further, note that (2b) is an indirect 
confirmation that the specific heat is linear in T as it is for the non-interacting 
Fermi gas. 



IV. Test of Fermi Statistics 

We have carried out Molecular Dynamics simulations to test the ther- 
modynamic theory of ref.[l] by specifically measuring the density profile, its 
center of mass and the fluctuations of the center of mass. The MD code 
used in this paper was provided by Jysoo Lee. It is identical to those used 
by HLRZ group in Germany [8] and thus we refer to ref. [8] for details about 
the code. To start a simulation, the program first places the particles into 
a two-dimensional configuration specified in the code by the user. Then 
the particles are allowed to relax under gravity and the particle-particle and 
particle- wall interactions before any vibration is "turned on" . The vibration 
of the bottom plate should not begin until the particles are motionless. This 
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relaxation time is found through trial and error while tracking the motion of 
the particles. The configurational average was taken every 100 time steps, 
and t N 2 • 10 7 , which is approximately 2000 vibrating cycles. Then, the 
center of mass and the fluctuations are plotted as a function of time. We 
now present the data. 

1. The Fermi Temperature T 

Simulations were carried out in two dimensional boxes with vertical walls 
for N particles each with a diameter D = 0.2cm and a mass m = 4n(D/2) 3 /3 = 
4.188»10- 3 #m with degeneracy Q, using (N,fi) = (100,4), (200,4), and (200,8) 
with a sine wave vibration and (100,4) and (200,4) with a triangular wave. 
Here, degeneracy is simply the number of columns. The dimensionless Fermi 
energy /j, — \ij D — N/Q, which is the initial layer thinckness. Intially, the 
particles were arranged on a square lattice, and the frequency cu was set to a 
fixed value, uo = 2nf = AOtc, which is the experimental value (f=20Hz) used 
by Clement and Rajchenbach [3]. The strength of the vibration, T = Auj 2 /g, 
was varied by changing the amplitude, A, of the vibration. Initially, T was 
varied from 1.2 to 20 in the simulations with g = 981cm/ s 2 . The density 
profile closely follows the Fermi profile for low T's, but it crosses over to 
the Botzmann distribution for higher T's (Fig.l). Note that the deviation 
from the Fermi profile for large T is quite noticeable. The deviation of the 
T = 20 data from the Boltzmann function near the bottom plate is due to 
the model used to approximate particle-particle reactions in the simulation 
program. The springs used in the soft sphere model of these particles cause 
the odd motion observed in the bottom layers of the columns of particles. 
These bottom particles have more of a spring force than the top layers be- 
cause the extra mass on top causes the bottom particles to deform more. 
This deformation leads to the separation of these bottom particles from the 
bulk, giving the density profile its odd shape from < z < 4 for T = 20. 
Hence, our analysis is focused strictly for T < 4. A/D changes from 0.373 
for T = 1.2 to 1.24 for T = 4. For all the data, 4.90 < R = fxD/TA < 111.20, 
and thus, the Fermi statistics is satisfied (See ref. [1] for this criteria). 

Representative samples of the density profiles for the sine wave with 
(N,f2) = (100,4) for T = 1.2,1.5,2.0,2.5,3.0,3.5 and 4.0 are plotted in Fig.2. 
Similar profiles were obtained for other sets of parameters. The data are then 
fitted by the Fermi function, and the Fermi temperature was determined. 
(Note that the temperature, Tf, discussed in the following text is actually 
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T/mg. This is because the fitting function p(z) = p c /(l + exp[(z — p)/Tf}) 
(See eq.(2)) was used with p, = p Q D and z — (i — 1/2)D (height with units 
of centimeters.) 

If h (T) with a = 64/5 is used in our temperature formula, the predicted 
values match extremely well with the measured ones. A reasonable agreement 
between the theory and the simulation is reached when one uses h (T) for 
both the sine and triangular wave data as shown in Fig.3a and 3b. This 
agreement is fairly constant for all five simulated systems. Now we turn our 
attention to the study of the temperature scaling. 

2. Study of temperature scaling 

We now test the scaling relation predicted by eq.(la) between the tem- 
perature and control parameters such as frequency (u), gravitational accel- 
eration (g), and particle diameter (D). 

(a) Frequency Dependence: The frequency of the sine wave was varied 
for fixed values of T = 2.0, N=200, and Q = 4, over the range 5 < uo < 35. 
The temperature was found to scale with frequency as 

where m 1.16, which is fairly close to the predicted value mi — 1. (Fig. 4a) 

(b) Gravity Dependence: Temperature as a function of the gravita- 
tional acceleration is displayed in Fig.4b with fixed values of T = 2.0, N=200, 
Q, = 4, and / = u/2n = 20Hz for the range of 150 < g < 1200. The tem- 
perature scales as 

T « g~ m2 

where m 2 ~ 0.48, which is again, close to the predicted value m 2 = 0.5. 
We notice a change in the temperature dependence for g c > g with g c m 
1000cm/s 2 , the reason for which is not so clear at this point. 

(c) Diameter Dependence: The temperature as a function of the di- 
ameter, D, at the fixed frequency / = uo/2-k = 20Hz and g = 981cm/ sec 2 is 
displayed in Fig. 4c with fixed values of T = 2.0, N=200, and Q = 4 for the 
range 0.025 < D < 1.60. Here, as the diameter was changed, the density was 
changed accordingly using p = (3/A7rD 3 )m to ensure that the mass remained 
at a constant value of 4.188a;10" 3 gm. The temperature scales as 

T w D~ m ' 6 
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where m^, m 0.53, which is again, close to the predicted value for = 0.5. 

To summarize, we find that the scaling relations between the temperature 
and u, g, and D appear to be consistent with the theoretical predictions 
within the error bars. Our next task is to examine the scaling relations 
presented for the center of mass and its fluctuations. 

3. Center of Mass 

The center of mass is plotted in Fig. 5 as a function of T 2 in a graph 
of seven different T's for N=200, Q = 4 and Q = 8, and N=100, Q = 4 
under a sine wave vibration and N=100, Q = 4 and N=200, Q = 4 under the 
triangle wave vibration. The solid line, a guide for the eye, seems to confirm 
the scaling predictions of the Fermi statistics as given by the formula (2a), 
namely 

Az oc T 2 . 

However, there exists a large discrepancy in the amplitude, C. For N=200 
and fl = 4, the theory predicts that one should get from (2a) C = n 2 /6fi = 
7T 2 /6// ~ .16 with jj, = Dfi Q , the actual height of the Fermi surface. The 
simulation results yield, C ~ 3.0 for the sine wave and C ~ 5.6 for the 
triangular wave, a difference of a factor 20 and greater when compared to 
with the theoretical prediction. This discrepancy has been traced in some 
detail using Mathematica, and it appears that the center of mass defined in 
(2a) is extremely sensitive to a slight change in the Fermi energy, /x. For 
example, if the average fi = 10.01 is used for the N=200, Q = 4 sine wave 
data and the integral is computed with Mathematica, the plot of the center 
of mass values as a function of T 2 yield a slope that is indeed 0.16, just as 
predicted by the Fermi analysis. However, if the experimentally determined 
ranging from 9.91 to 10.13 is used and the center of mass is computed using 
Mathematica, the slope turns out to be 2.5. Hence, one can conclude that 
the temperature dependence of /i causes the huge discrepancy observed in the 
amplitude. Note that when the density of states is independent of the energy, 
the Fermi energy must remain constant. Yet, in real simulations, it changes 
for each value of T. Fortunately, this 5 % change in the Fermi temperature 
does not make a noticeable difference in the Fermi fitting. Therefore, the 
temperature fitting is robust. 

4. Fluctuations of the Center of Mass 
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The fluctuations of the center of mass were also measured and plotted as 
a function of T 3 for N=200, 0, = 4 and fl = 8 and N=100, ft = 4 under a 
sine wave vibration and N=200, = 4 and N=100, SI = 4 under a triangle 
wave vibration. (Fig. 6) Once again, this confirms the validity of the Fermi 
statistics, which implies that 

< (Az) 2 >oc T 3 . 

The proportionality constant, though, is a way off from the theoretical value, 
Dir 2 /3fi 2 ~ 6.4 • 10~ 3 . The actual slope is of the order 1. Considering the 
fact that the amplitude in the center of mass is off by a factor 10, one can 
expect the errors in the fluctuations to be larger, of the order 10 2 . This is 
again due to the sensitivity of the Fermi integral to /i. A small change in 
H is magnified in the fluctuations, < {Az 2 ) >, resulting in the 10 2 factor 
difference. An additional source of discrepancy is due to the fact that the 
fluctuations of the center of mass are quite large in the vibrating bed, and 
all the particles down to the ones at the very bottom of the layer fluctuate 
in a continuum space. This is in contrast to the Fermi particles of a lattice 
model, where most of the particles below the Fermi surface are locked and 
are thus, inactive. Hence, while the average position of the grains appear to 
obey the Fermi distribution function quite well, its actual magnitude in the 
fluctuations may not. The surprise here, however, is that the temperature 
dependence of the flutuations still appears to obey the T 3 law of the Fermi 
statistics. 

V. Discussion 

There are several issues that should be discussed in connection with the 
present work. 

First, the thermodynamic aspect of the theory presented in [1] should be 
investigated in more detail in the future. Certainly the temperature formula, 
eq.(la) and (lb), certainly needs to be modified and improved. Further, note 
the obvious fact that the temperature in the Fermi distribution function is 
not the same as the kinetic temperature, T k = m < v 2 > /2. T k is zero 
when there is no motion, but the temperature T defined in ref [1] is non zero 
even though there is no motion because of the entropic contribution. In this 
sense, our theory is similar to that proposed by Edwards [9]. Edwards and 
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his collaborators proposed a thermodynamic model for an equilibrium state 
by defining a temperature like variable, termed the compactivity, 




where V is the total volume of the granular system, and S is its entropy. 
The volume V can be determined easily, but there is no systematic way 
of determining the entropy, S, of a disordered system. Hence, there is no 
calculational mechanism to relate the central parameter of his theory, the 
compactivity X, to the experimentally controlled parameters such as T. Thus, 
in our opinion the theory is not closed. Within his formalism, it is quite 
a challenge to systematically compute the entropy of a disordered system. 
Recently, the Chicago group [10] determined X using a fluctuation formula 
such as eq.(2b). Note, however, that this is an experimental method used to 
determine X, and a formula relating X to control parameters is still missing. 

Our theory departs from Edwards in two ways: First, it is a thermody- 
namic theory of the configurational statistics of a nonequilibrium dynamical 
state or steady state. In this view point, the dynamics of compaction [10,11] 
might be viewed as a transition of an unstable state into an equilibrium 
state [12]. The temperature equation may be ignored in the hydrodynamic 
approach of the convective instability of granular materials [13,14]. Next, 
the relationship between the temperature variable and the control parame- 
ters was determined by comparing the kinetics and thermodynamics. This 
provides a closure to Edward's thermodynamic theory by providing a spe- 
cific relationship between the temperature or compactivity and the external 
control parameters. 

Second, the thermodynamic theory in the presence of gravity with a hard 
sphere gas presents some interesting puzzles and needs to be pursued in 
the future, notably the question of whether or not the temperature is an 
extensive quantity along the direction of anisotropy. For example, consider 
point particles confined in a three dimensional box of size (L x , L y ,L z ) under 
gravity along the z direction, for which the energy level is given by e = 
p 2 /2m + mgz. One can show easily that the energy per particle u is given by 



u = E/N 



1 dlnZ 
N dp 



\kT - f{x) 



(3) 
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where Z = z N /N\ is the partition function of the N particle system, z = 
(2mir/P) 3 / 2 L y L x (l—exp(—(3x))/f3mg, &ndf(x) = xexp(— (3x)/(l— exp{ — (5x)) 
with x = mgL z . In the limit L z — > oo, the energy per particle obeys the 
equipartition theorem and approaches hkT/2. Further it approaches zero at 
the zero temperature limit. Note also that the total energy E is an extensive 
quantity, and the thermodynamic relation between the temperature and the 
entropy also is satisfied, namely: S/k = InZ + j3E and 



However, for granular materials with finite diameter D(hard spheres), the 
situation becomes a little different. In the zero temperature limit, while the 
kinetic energy approaches zero, the potential energy is still a function of dis- 
order and it is not zero. For point particles, both the kinetic and potential 
energy approach zero because particles can be compressed indefinitely. Fur- 
ther, for a hard sphere gas, the potential energy is not extensive and thus the 
temperature T is not intensive. Consider for example the granular materials 
confined in a two dimensional box of width L and the height H, for which the 
potential energy E = LHf = LH 2 /2. Now, change the box size by a factor 
two: if the width is increased to 2L, then the potential energy doubles, i.e: 



Hence, the energy is extensive. On the other hand, if the height is increased 
by a factor 2, then the potential energy increases by a factor 4, i.e: 



So, the energy is not extensive. Hence, one finds here that while the energy 
does not depend on the way the system is doubled for point particles, it 
does for hard spheres with finite diameter. The thermodynamic definition 
of the temperature with hard spheres in the presence of gravity does not 
seem so simple as in the case of point particles. This is not inconsistent with 
recent experimental result [15], where it was reported that even though the 
velocity profiles obey perfect Gaussian, one needs two different temperatures 
to describe the velocity profile of excited grains in a vibrating bed 





E = 2LH • — = 2E 



E = L»2H* = AE . 
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Third, comments can be made about the MD code which was used. It 
is important to note that the MD code used in this work allows the grains 
to deform upon contact and the amount of deformation, in particular along 
the normal direction, depends on the mass and the spring constant. Within 
this MD code, unlike the assumption made in the theory that particles are 
compact and nondeformable, the temperature does depend weakly on the 
mass of the grain. More precisely, when a spring with the spring constant 
K is displaced a distance Az by a mass m, the velocity of the mass is given 
by \mv 2 = ^K(Az) 2 , i.e. v = (K/m) 1 / 2 . Now, note that in the MD code 
[8] used in this paper, the normal contact force at the bottom plate, F, has 
three components. 

F = K n (Ar) m -mg- ~f n m eff (v • n) 

where the first term is the Hertzian contact force [16], the second term is 
the gravity term, and the third term is due to the dynamic friction. Since 
the acceleration of the particle, a, is given by a = dv jdt = F n /m, and the 
first term is independent of the mass, the bouncing velocity of the particle 
at the bottom plate is a function of the mass for a given K n . Therefore, the 
jumping height of the ball does depend on the mass of the particle. While 
the MD codes that allow such deformation are based on certain models [17] 
and experiments [18], it is not certain at this point whether such mass depen- 
dent dynamics of deformable grains is a realistic modeling of real granular 
materials. Presumably, this ambiguity might have had some effect on the 
discrepancies in the amplitude for the center of mass and its fluctuations. 
It is highly desirable to carry out the same study undertaken here with a 
hard sphere MD code [19] which does not allow such deformations to occur. 
Finally, we also point out that the microscopic basis of the Fermi statistics 
has been recently examined by one of the authors in ref. [20], where it was 
demonstrated that the crossover from Boltzmann to Fermi statistics as the vi- 
brational strength decreases arises via the condensation of granular particles 
under gravity. This is another consequence of the excluded volume inter- 
action of finite grains that one cannot compress such a system indefinitely, 
which is the key to uncovering an interesting phenomena of 'condensation of 
grains under gravity. Putting aside the mathematical aspect of this conden- 
sation, which is mathematically similar to Bose condensation, the underlying 
physics of this phenomenon is not difficult to understand. If we start with 
particles at a high temperature and decrease the temperature slowly, then a 
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portion of particles begins to settle down from the bottom, and two phases 
develop: a fludized region near the top and a solid or glass region near the 
bottom. This is what we term the condensation phenomenon of granular par- 
ticles, and in the presence of such two phases, the continuum theory may pose 
some problems. For detailed summary of how the kinetic theory describes 
such a condensation phenomenon, see ref.[20]. 
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Figure Captions 



Fig.l. Density profile for T = 20 (square) and T = 2 (circle) with N=100 and 
Q = 4. For T = 2, the inequality eq.(3) that insures the validity of the Fermi 
statistics is satisfied. The Fermi statistics breaks down for T = 20 and the 
density profile is approaches that of the dilute gas with gravity, shown with 
the dashed line. 

Fig. 2 Representative density profiles for different T 's are plotted for N=100, 
tt = 4. 

Fig. 3. Comparison between the measured temperatures by MD and the pre- 
dicted ones by eq.(2a)and (2b) for (a) sine wave excitations (b) triangular 
wave excitations. 

Fig. 4. The log-log plot of the scaling of the measured temperature, T, against 
the controlled parameters (a) frequency - co, (b) gravity - g, and (c) the 
diameter - D. The solid lines are a guide for the eye and their slopes are: 
1.167(w), 0.48(g) and 0.528(D) respectively. 

Fig.5. The center of mass, < z(T) >, is plotted for N=100 Q = 4 and N=200 
Q = 4 for both the sine wave and triangle wave vibration and N=200 Q = 8 
for the sine wave vibration as a function of T 2 . The straight lines are a guide 
to the eye. 

Fig. 6. Fluctuations in the center of mass, < (Az(T) 2 >,is plotted for N=100 
Q = 4 and N=200 Q = 4 for both the sine wave and triangle wave vibration 
and N=200 VI = 8 for the sine wave vibration as a function of T 3 . The 
straight line is a guide to the eye. 
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